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Abstract 

Recent discoveries of direct acting antivirals against Hepatitis C virus (HCV) have raised hopes of 
effective treatment via combination therapies. Yet rapid evolution and high diversity of HCV 
populations, combined with the reality of suboptimal treatment adherence, make drug 
resistance a clinical and public health concern. We develop a general model incorporating viral 
dynamics and pharmacokinetics/pharmacodynamics to assess how suboptimal adherence affects 
resistance development and clinical outcomes. We derive design principles and adaptive 
treatment strategies, identifying a high-risk period when missing doses is particularly risky for de 
novo resistance, and quantifying the number of additional doses needed to compensate when 
doses are missed. Using data from large-scale resistance assays, we demonstrate that the risk 
of resistance can be reduced substantially by applying these principles to a combination 
therapy of daclatasvir and asunaprevir. By providing a mechanistic framework to link patient 
characteristics to the risk of resistance, these findings show the potential of rational treatment 
design. 



Introduction 

Hepatitis C virus (HCV) affects approximately 170 million people and chronic infections can lead 
to cirrhosis and hepatocellular carcinoma (Lavanchy 2009; Thomas 2013). Recently, development 
of direct acting antivirals (DAAs) against HCV infection has revolutionized the field of HCV 
treatment, because of their high potency, broad applicability and mild side effects (Gane 2011; 
Scheel and Rice 2013). Combination therapies of DAAs have achieved remarkably high rates of 
sustained virological response in clinical trials (Lok, Gardiner et al. 2012; Pol, Ghalib et al. 2012; 
Afdhal, Zeuzem et al. 2014; Feld, Kowdley et al. 2014; Kowdley, Gordon et al. 2014; Zeuzem, 
Jacobson et al. 2014). However, due to the relatively low genetic barriers of most DAAs 
(Pawlotsky 2011; Robinson, Tian et al. 2011; Qi, Olson et al. 2014), the high intrinsic mutation 
rate of HCV (Powdrill, Tchesnokov et al. 2011; Ribeiro, Li et al. 2012), and the high viral diversity 
(Pybus, Charleston et al. 2001; Simmonds 2004; Thomas 2013), combined with the reality of 
suboptimal treatment adherence (Lo Re, Teal et al. 2011; Lieveld, van Vlerken et al. 2013), viral 
resistance represents a clinical and public health concern (Sarrazin and Zeuzem 2010; Pawlotsky 
2011). Indeed, resistance to single DAAs has already been observed frequently for many 
candidate DAAs, and patients must be treated with combination therapies to avoid treatment 
failure. If not properly managed, resistance could quickly develop to combination therapies and 
render these new DAAs useless, as observed for other antimicrobial treatments, squandering the 
potential health gains from these recent breakthroughs (DiMasi, Hansen et al. 2003; Roberts, 
Hota et al. 2009; Smith, Okano et al. 2010). 

Suboptimal patient adherence to dosing regimens is a crucial risk factor for resistance 
development in both HIV and HCV treatments (Paterson, Swindells et al. 2000; Bangsberg, Perry 
et al. 2001; Lo Re, Teal et al. 2011; Lieveld, van Vlerken et al. 2013). Although high rates of 
sustained virological response have been achieved in clinical trials, adherence levels may vary 
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substantially among the vast population of infected patients, owing to long treatment periods 
and complicated regimens associated with DAA combination therapies (Weiss, Brau et al. 2009; 
Lo Re, Teal et al. 2011; Evon, Esserman et al. 2013; Gordon, Yoshida et al. 2013; Lieveld, van 
Vlerken et al. 2013). Rational design of combination therapy regimens, enabling individualized 
regimens based on the genetic composition of a patient's infection and real-time adjustments for 
missed doses, is a top research priority to avoid resistance (Spellberg, Guidos et al. 2008; Gelman 
and Glenn 2010; zur Wiesch, Kouyos et al. 2011; Lieveld, van Vlerken et al. 2013). Mathematical 
models are well suited to address this problem. Previous modeling studies for HIV infections 
have illuminated potential mechanisms underlying treatment failure and explained puzzling 
clinical observations (Wahl and Nowak 2000; Rosenbloom, Hill et al. 2012). However, HCV is a 
curable disease and its infection, goal of treatment and mechanism of resistance differ from HIV 
in many respects (Soriano, Perelson et al. 2008), including no known latent reservoir and a finite 
treatment period to eradicate the virus. Here, by integrating 
pharmacokinetics/pharmacodynamics (PK/PD) and viral dynamics into mathematical models, we 
develop the first general theory to assess the impacts of suboptimal adherence on the outcome 
of DAA-based therapies for HCV infection. We derive design principles that can be generalized to 
therapies involving different classes and different numbers of drugs. Using large-scale data from 
in vitro resistance assays and human clinical trials, we apply this framework to a combination 
therapy of daclatasvir and asunaprevir (Suzuki, Ikeda et al. 2012), and derive evidence-based 
adaptive treatment strategies for treatment protocols over time according to resistance profiles 
and adherence patterns. 



Results 



Resistance to antiviral treatments can develop through selection of preexisting mutants or de 
novo generation of new mutants. A core principle for designing effective combination therapy is 
that, if patients fully adhere to the treatment regimen, the treatment must suppress all 
preexisting mutants and de novo resistance should be unlikely (Ribeiro and Bonhoeffer 2000). 
Missing doses, however, can lead to suboptimal drug concentrations, allowing growth of some 
preexisting mutants with partially resistant phenotypes. Growth of these mutants allows the viral 
population to survive longer, possibly generating further mutations that contribute de novo 
resistance against the full combination therapy. For example, consider a combination therapy of 
two DAAs, A and B, as shown in Fig. 1A. If missed doses and pharmacokinetics lead to a drop in 
the concentration of drug A, this allows growth of the preexisting mutant, AB', (which is already 
resistant to drug B), thus opening opportunities to generate the fully resistant mutant, A'B'. 
Therefore, the dynamics of the subset of preexisting mutants that have a high level of resistance 
against single DAAs determine resistance evolution and treatment outcomes for combination 
therapies. In the following, we denote these mutants as 'partially resistant' mutants. 



The effective viral fitness, R e ff(t) 

The fitness of a particular strain in a treated patient is determined by the PK/PD of the drug, the 
level of resistance of the strain, and the availability of target cells, i.e. uninfected hepatocytes for 
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HCV (Fig. IB). We can integrate all these factors (for any class of DAA therapy) into a single 
number, the effective reproductive number under treatment, R e /f{t) (Fig. 1C). R e ff(t) is the average 
number of cells infected by viruses produced by a single infected cell. It acts as a measure of viral 
fitness, and can be calculated as: 

R eff (t) = (\-e(r))-R 0 -h(t) (l) 

where t is time since treatment starts, r is the time since last dose, e(t) is the efficacy of the drug 
combination at time x during the dosing cycle, R 0 is the reproductive number of the virus in the 
absence of treatment, and h{t) is the normalized abundance of target cells (see Supplementary 
Materials). Under effective treatment, the availability of target cells, h{t), increases quickly to 
reach the infection-free level (Rong, Dahari et al. 2010), and therefore, the overall viral fitness 
increases over time as h{t) increases under effective treatment (Fig. 1B,C). When adherence is 
optimal, the value of R e ff for a partially resistant mutant is always less than 1 (i.e. viral 
suppression); however, if doses are missed, drug concentration declines exponentially and R e ff 
can become greater than 1 (i.e. viral growth) (Fig. 1C). 

The growth of partially resistant mutants and the need for extended treatment 

We now consider how suboptimal adherence impacts the dynamics of partially resistant mutants. 
As an illustration, we contrast simulations assuming perfect adherence versus suboptimal 
adherence. Missing doses leads to rapid decreases in drug concentration, and thus increases in 
R e ff of a partially resistant mutant (Fig. 2A-C). This means that extra doses are needed to 
compensate for the missed doses to suppress the mutant to extinction (Fig. 2D), and also that 
the number of newly infected cells rises substantially, which increases the opportunity for de 
novo resistance (Fig. 2E). 

We approximate the time-varying values of R e ff{t) during periods when doses are missed, by 
calculating the average effective reproductive number, R aV e,mj as (see Materials and Methods): 

K V e, m (t)~Q-£ ave , m )-R 0 -Kt) (2) 

where t is the time when the patient starts to miss doses, m is the number of consecutive doses 
missed and e me , m is the average drug inhibition during the period when m consecutive doses are 
missed. This allows us to generalize our theory to any DAA combinations for which £ ave , m can be 
either estimated from pharmacokinetics/pharmacodynamics data or calculated from mutant 
resistance profiles (Wahl and Nowak 2000). 

We then ask, if m consecutive doses are missed beginning at time t, how many extra doses, N m , 
are needed to compensate? This number, which we denote 'compensatory doses', can be 
approximated as (see Materials and Methods): 

N m (t)~m- R ave>m (t) ~ m • (1 - e me>m )-R 0 - hit) (3 ) 

This allows us to estimate the total duration of treatment needed to clear infection for a given 
adherence pattern. Furthermore, since h{t) increases over time under effective treatment (Rong, 
Dahari et al. 2010), Eqn. 3 shows that a higher number of extra doses are needed to eliminate 
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the infection if doses are missed later in treatment. 
De novo generation of fully resistant mutants 

To assess the risk that a partially resistant lineage will give rise to full resistance, we calculate the 
expected number of target cells, <£ m , that become infected by fully resistant mutant viruses due 
to de novo mutation during a period when m consecutive doses are missed. This quantity is the 
product of the cumulative number of cells newly infected by a partially resistant mutant and the 
effective mutation rate from that mutant to the fully resistant mutant, //^(see Materials and 
Methods): 

6 (Q ^^^^ 

(0 - H« ■ Of) ■ - ■ (e^^ T - 1) (4) 

cumulative number of cells newly infected by a partially resistant mutant 

where /(f) is the number of cells infected by the partially resistant mutant at time f when the first 
dose is missed, and 0(f) represents the potential to generate new infections. 5 is the death rate 
of infected hepatocytes, and T is the scheduled interval between two doses. O m quantifies the 
risk that a fully resistant mutant infects target cells, but whether it emerges and becomes 
established within the host depends on its fitness and the stochastic dynamics of invasion 
(Alexander and Bonhoeffer 2012; Loverdo, Park et al. 2012; Loverdo and Lloyd-Smith 2013). 

The strong dependence of <5 m on //^predicts that designing combination therapies to increase 
the genetic barrier to full resistance, e.g. using DAAs with higher genetic barrier or adding an 
extra drug into the combination, can reduce <b m by orders of magnitude or more, thus it would 
lead to drastic reductions in the probability of generating full resistance (compare trajectories a 
and b in Fig. 3A). 

Eqn.4 also allows us to assess when during treatment it is most risky to miss doses, which can 
inform treatment guidelines. Changes in two quantities, /(f) and 0(f), determine changes in O m 
over the course of a treatment regimen. For as long as adherence is perfect, /(f) decreases 
exponentially, while 0(f) increases over time since R me ,m[t) increases as the abundance of target 
cells rises over time (Fig. 3B). Thus the value of 4> m first increases (due to rapid recovery of target 
cells) and then decreases exponentially (due to decrease of infected cells). This leads to a 
high-risk window period, during which missing doses is especially risky for generating full 
resistance (Fig. 3A). This qualitative finding is robust to changes in model parameters, though 
quantitative predictions of the risk of full resistance depend on the fitness of the mutant (/? 0 ), the 
half-life of infected cells (8), and the rate at which the target hepatocytes become available upon 
treatment (Fig. SI). 

Design principles and adaptive treatment strategy for DAA combination therapy 

These results suggest principles for rational optimization of treatment outcomes. Individualized 
therapies could be designed for patients with risk factors for low adherence, by selecting drug 
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combinations that impose a higher genetic barrier than required to suppress all preexisting 
mutants, to reduce the risk of de novo resistance. 

Adaptive treatment strategies could be developed based on the theoretical findings shown 
above. For a particular combination therapy, the high-risk window period for missing doses can 
be calculated by integrating the values of <&,„ for all partially resistant mutants present in a 
patient. Then, for patients with risk of low adherence, supervised dosing during the high-risk 
window period would reduce the risk of resistance and treatment failure. Another alternative is 
to treat the patient using a higher number of DAAs in combination during the high-risk period, 
and then switch back to a combination with a lower number of DAAs afterwards. If doses are 
missed during treatment, the patient should be treated with extra doses, computed as the 
maximum value of the N m values calculated for all partially resistant mutants. For the lowest risk 
of de novo resistance, the prescribed number of compensatory doses {N m ) should be taken, 
uninterrupted, immediately after doses are missed. Otherwise the infected cell population may 
rebound to a high level, which can make further missed doses very risky for resistance. 

Case study: combination therapy ofdaclatasvir and asunaprevir 

To demonstrate the practical applicability of our theory, we consider a recently developed 
interferon-free combination therapy based on an NS5A inhibitor, daclatasvir, and an NS3 
protease inhibitor, asunaprevir (Suzuki, Ikeda et al. 2012). In clinical trials, a large proportion of 
patients infected with HCV genotype-lb achieved sustained virological response (i.e. viral 
eradication) when treated with daclatasvir and asunaprevir for 24 weeks, although viral 
breakthrough and viral relapse occurred in a small fraction of patients (Karino, Toyota et al. 2013; 
Kosaka, Imamura et al. 2014). 

We first consider patients with the wild-type virus at baseline, i.e. the wild-type virus is the 
dominant strain before treatment. Using the PK/PD data for each drug (Eley, Pasquinelli et al. 
2010; Nettles, Gao et al. 2011; Ke, Loverdo et al. 2014) and the resistance profiles data measured 
for genotype-lb HCV (Fridell, OJu et al. 2010; Fridell, Wang et al. 2011), we predicted which 
mutants are potentially fully-resistant to this combination therapy and calculated the values of 
N m and <!>,„ for each of the partially resistant mutants (Fig. 4A,B) (see Supplementary Materials 
for more detail). Choosing the highest values of N m and <£>,„ among all the partially resistant 
mutants allows us to project the overall risk arising from missed doses over the course of 
treatment, and we found required numbers of compensatory doses were modest and the risk of 
de novo resistance is low (Fig. S2A). To demonstrate that the theoretical approximations 
represent the full viral dynamics accurately, we simulated a multi-strain viral dynamics model 
(see Materials and Methods), assuming 1-3 day blocks of consecutive doses are missed randomly 
within a treatment regimen lasting 24 weeks. The model predicts that relapse of L31M+Y93H or 
L31W would be observed when overall adherence is less than 90% (Fig. 4C,D). Indeed, the 
L31M+Y93H mutant has already been detected in one relapse patient in a clinical trial (Karino, 
Toyota et al. 2013). There is excellent agreement between simulation results and theoretical 
predictions (based on Eqn.3 and 4) for the number of cells infected by different mutants after 24 
weeks of treatment and the cumulative number of cells infected by partially resistant mutants 
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over the treatment period (Fig. 4D and S3). 

We then simulated outcomes when the doses are guided by the adaptive treatment strategy 
{guided dosing). Because the risk of de novo resistance when doses are missed is low, there is no 
high-risk period for de novo resistance in this case (Fig. 4B). If patient dosing is guided, i.e. all the 
required doses and the extra doses to compensate for the missed doses are taken, the infection 
can be cleared successfully (Fig. 4E). Again, we find excellent agreement between simulation 
results and theoretical predictions (Fig. 4F). 

Many patients bear the Y93H mutation at baseline and this mutation reduces the genetic barrier 
to full resistance by one nucleotide(Karino, Toyota et al. 2013). Our theory suggests that 
reducing the genetic barrier to full resistance will drastically increase the risk of treatment failure. 
We repeated our analysis for patients with Y93H at baseline, to test how our adaptive treatment 
strategy works when the risk of resistance is high. As predicted, many more days of treatment 
are needed to compensate for missed doses, and the risks of generating full resistance de novo 
are high (>0.01) during the first 3 weeks of effective treatment if 2 consecutive doses are missed 
(or first 4 weeks if 3 doses are missed; Fig. 5A,B and S2B). De novo full resistance is likely if doses 
are missed randomly and adherence is less than 90% (dark red area in Fig. 5C). The predicted 
number of infected cells agrees well with simulation, except when adherence is very low such 
that viral load rebounds back close to the pre-treatment level (Fig. 5D and S4-S6). In stark 
contrast, when doses are guided, the risk of de novo resistance becomes much lower (compare 
Fig. 5C with 5E). Again, for patients who do not clear infection after 24-week treatment, 
extended periods of treatment as predicted by our theory (using Eqn.3) can clear infection with 
low risk of resistance. The efficacy of the adaptive treatment strategy is robust across different 
parameter values (Fig. S7-S12 and Supplementary Materials). Therefore, our treatment strategy 
can improve clinical outcomes substantially by adjusting on-going treatment based on patient 
adherence patterns. 



Discussion 

With the prospect of interferon-free combination therapies becoming available to the HCV 
infected population (Lok, Gardiner et al. 2012; Pol, Ghalib et al. 2012; Scheel and Rice 2013; 
Thomas 2013; Afdhal, Zeuzem et al. 2014; Feld, Kowdley et al. 2014; Kowdley, Gordon et al. 2014; 
Zeuzem, Jacobson et al. 2014), there is an urgent need to design treatment strategies that will 
prevent or delay the development of resistance to DAAs. Extensive laboratory efforts have 
characterized the PK/PD parameters and mutant resistance profiles of DAAs (Eley, Pasquinelli et 
al. 2010; Fridell, Wang et al. 2011; Nettles, Gao et al. 2011; McPhee, Friborg et al. 2012; Oj, Olson 
et al. 2014). In this study, we integrate PK/PD parameters and viral dynamics into a unified 
framework to assess the impacts of suboptimal treatment adherence on the risk of treatment 
failure. This framework also enables adaptive management of DAA treatments. Using simulations 
incorporating PK/PD and resistance profile data collected previously (Fridell, Qiu et al. 2010; 
Fridell, Wang et al. 2011; Nettles, Gao et al. 2011), we showed that treatment outcomes of 
combinations therapies of daclatasvir and asunaprevir can be greatly improved by this adaptive 
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treatment strategy, especially when the Y93H mutant is the dominant strain before treatment 
begins. 

For therapies with low genetic barriers to resistance, we have identified a high-risk window 
period during which de novo resistance is likely if doses are missed. Intervention efforts should 
focus on enhancing patients' adherence during this period. Additional complementary strategies 
could further reduce the risk of treatment failure. First, if doses are missed during the high-risk 
window, the immediate addition of another drug with a different mechanism of action from 
existing drugs may eliminate any low level of fully resistant mutants that has arisen. Alternatively, 
a patient could be treated preemptively using additional drugs during the entire high-risk period 
and switched to fewer drugs afterwards. Our theory also predicts the number of compensatory 
doses {N m ) needed to compensate for missed doses, in order to eliminate preexisting mutants. 
Interestingly, clinical trials have shown that adherence levels tend to decrease over time (Weiss, 
Brau et al. 2009; Lo Re, Teal et al. 2011); we show that more doses are needed to compensate 
for missed doses that occur later in treatment because of the rebound of target cells. Overall, 
these results highlight the importance of viral genotype screening and adherence monitoring. 
While many previous studies have focused on average adherence (Wahl and Nowak 2000; Weiss, 
Brau et al. 2009; Lo Re, Teal et al. 2011; Evon, Esserman et al. 2013; Gordon, Yoshida et al. 2013; 
Lieveld, van Vlerken et al. 2013), we emphasize that the timing of the missed doses is also a 
critical determinant of treatment outcome and the risk of resistance. 

There exist substantial heterogeneities among patients owing to variation in HCV genotypes, 
patient viral loads, death rates of infected cells (Neumann, Lam et al. 1998; Rong, Dahari et al. 
2010) and effectiveness of drug penetration (Ke, Loverdo et al. 2014). Our analysis has identified 
several factors that influence the impact of suboptimal adherence, particularly the rebound rate 
of target cells under treatment, the half-life of infected cells and the overall viral fitness, R 0 . We 
used the best available estimates of these parameters, but further empirical work is needed. If 
resistance profiles and viral parameters could be measured directly from a specific patient, then 
our framework linking these factors could be tailored to give patient-specific guidelines. 

Certain model assumptions reflect uncertainties in our current knowledge of HCV infection. First, 
our prediction about time to viral extinction should be treated cautiously. We predict the time of 
extinction (as in other models (Snoeck, Chanu et al. 2010; Guedj and Perelson 2011)) by 
assuming that infected cells decline at a rate set by their death rate, and infection is cleared 
when the number of infected cells is below one. However, factors such as pressures from the 
immune system and infections in different tissue compartments may influence the extinction 
threshold. Furthermore, if DAA treatment causes intracellular viral RNA to decay with negligible 
replication (Guedj, Dahari et al. 2013), the decline of infected cells may result from a 
combination of cell recovery and death of infected cells. Indeed, sustained virological response 
has been observed in clinical trials of DAA combination therapies with shorter durations of 
treatment (Poordad, Lawitz et al. 2013). Our model can be adjusted easily once the decay 
dynamics of infected cells are understood better. Second, our model captures the main features 
of pharmacodynamics and viral dynamics by assuming quasi-equilibrium for viral populations and 
drug penetration into liver cells. Further work that incorporates detailed intracellular interactions 
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(Guedj, Dahari et al. 2013) and different body compartments may improve model accuracy, once 
pertinent parameters are measured. However, a more detailed model may become analytically 
intractable. 

This quantitative framework is a step towards developing a tool (for example, see Ref. (Garg, 
Adhikari et al. 2005)) for clinicians to design combination therapies and adaptively manage 
treatment regimens to achieve favorable clinical outcomes. It highlights the importance of 
characterizing resistance profiles of HCV, screening for resistant mutations before treatment, 
and monitoring adherence patterns during treatment, so that treatment can be designed and 
adjusted in an evidence-based manner. This framework can be adapted easily to combination 
therapies based on other DAA candidates, or treatments of other curable diseases without a 
reservoir. 



Materials and Methods 

HCV model and Viral fitness in the presence of drug, R e /f(t) 

To analyze the dynamics of the virus, we constructed an ordinary differential equation (ODE) 
model to describe the long-term within-host dynamics of a single HCV strain under drug 
treatment, based on an established model developed by Neumann et o/.(Neumann, Lam et al. 
1998) (see Supplementary Material). In the model, e represents the proportion by which the 
therapy reduces viral growth {s is in the range of 0 and 1). Then, the fitness of the virus, R e ff{t), is 
the product of the complement of the therapy's efficacy (l-e(t)), the reproductive number of the 
virus, R 0 , and the availability of target cells, h{t) (Eqn. 1). 

Average effective viral fitness when m doses are missed, R aV e,m 

To approximate the time-varying viral fitness, /? e //(t), during the period when m consecutive 
doses are missed, we assume that the abundance of target cells stays constant. This is a good 
approximation, because the length of the period when consecutive doses are missed tends to be 
short compared to the time scale of target cell rebound. Then the only time-varying quantity in 
Eqn. 1 is s(t). We can calculate the average level of drug inhibition during the period when m 
doses are missed, £ m e,m, by incorporating parameters for pharmacokinetics and 
pharmacodynamics (for example, see Wahl and Nowak(Wahl and Nowak 2000)). Then the 
time-average effective reproductive number, R m e,m(t), for a mutant when m consecutive doses 
are missed starting at time t can be expressed as Eqn. 2. In practice, because the precise number 
of target cells at time t is hard to estimate, we can approximate R me ,m by setting h{t)=l, and then 

Rave,m becomes R ave m (t) = (1 - E ave m ) • R 0 . Because h(t)<l, this always overestimates the viral 
fitness and thus is a conservative estimate in terms of guiding treatment. 

The number of compensatory doses needed (N m ) 

To calculate N m for each mutant, we make the simplifying assumption that the dynamics of the 
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viral populations are at quasi-equilibrium, because changes in the viral populations occur much 
faster than changes in infected hepatocytes. Then, the dynamics of the number of cells infected 
by mutant viruses, /(t), are described by: 

^p- = (R (t)-iys-i(t) (5) 

dt n 

where 5 is the death rate of infected hepatocytes. If we approximate R e ff{t) using the constant 
Rave,m for the period when doses are missed, Eqn. 5 can be solved analytically. Then, the number 
of infected cells after missing m consecutive doses starting at time f 0 can be expressed as: 

I(t 0 +m-T)~I(t a )- cxp((R avem (t 0 )- iyd-m-T) (6) 



We now consider the situation when m consecutive doses are missed, and ask how many 
uninterrupted doses (compensatory doses) must be taken so that the number of cells infected by 
the mutant is suppressed to a same number as if the m doses had not been missed. We first 
calculate the number of infected cells if the m consecutive doses are taken, i.e. if dosing is 
optimal: 

Iopumat (t 0 + m-T)~ I(t 0 ) ■ exp ((R ave>0 (f 0 ) - 1) • 8 ■ m • T) (7) 

where l{t 0 ) is the number of cells infected by the mutant at time r 0 , R 0 ve,o is the average effective 
reproductive number of the mutant when all doses are taken, and T is the scheduled interval 
between doses. 



We then analyze the situation where a patient skips m consecutive doses, starting at time r 0 , and 
then takes N m compensatory doses immediately afterwards. In this case, assuming the number 
of target cells does not change much during this period, we can approximate the number of cells 
infected by the mutant at the end of the N m doses as: 

(8) 

By equating the right hand sides of Eqn. 7 and 8 and solving the equation, we derive the 
expression for N m : 

NM- Rave - m(to) ~ Rave -° (to) -m (9) 
For potent therapies, usually R ave0 (t 0 ) = 0 . Then we get Eqn. 2. 



In the derivation above, we have assumed that the target cell abundance stays constant during 
the period under consideration. This would be a good approximation if only a few days of doses 
are missed or if the target cell has already rebounded to the infection-free level. If the 
abundance of target cells changes considerably during the period under consideration, an 
alternative, conservative approach would be to assume h{t)=l and take 

N mm3x (t 0 ) *=> m-(l-E ave m )- R 0 compensatory doses after missing m consecutive doses of 
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treatment. 

The number of doses to eradicate a mutant (N erad ) and the number of cells infected 
by a mutant (l(t)) 

One important application of N m is to predict the number of remaining doses needed to 
eradicate a mutant, N era d, in a patient during treatment. This number can be calculated as follows. 
If adherence is perfect, the number of infected cells declines exponentially at a rate set 

approximately by the death rate of infected cells, d: I(t) « I 0 • exp(-5 • t) , where / 0 is the 

number of cells infected by a mutant of interest before treatment. If we assume that a mutant 
goes extinct if the expected number of infected cells in a patient goes below 1, the number of 
doses needed to eradicate a mutant before treatment (assuming adherence is perfect), W era d,o, is 

calculated as: N m lo g( 7 o) . 

o- 1 

When doses are missed during treatment, it is clear from the calculation of N m above that N m -m 
extra doses of treatment are needed to eradicate the virus. Therefore, if a patient has taken a 
total of x doses and has had k instances of missing doses before time t, with m, days of doses 
missed in the / th instance (i=l,2,...,k), then the number of remaining doses needed to eradicate 
the mutant is calculated as: 

N erad =N eradfi -x + 2(N mJ -m i ) (10) 

i-l 

We can use Eqn. 10 to predict the number of cells infected by a mutant as: 
I(t) = exp((5- N erad (t)-T) . In our model, and a patient is cleared of infection when all mutants 

are driven to extinction. The accuracy of this approximation is shown in Figs. 4D,F and 5D,F. 



The risk of full resistance if doses are missed ( & m ) 

To calculate the risk of full resistance during the period when m doses are missed, we first 
calculate the number of cells newly infected by a partially resistant mutant when m doses are 
missed, fi m (f). Again, we use R ove ,m{t) to approximate R e fj{t), the total number of cells infected by 
the mutant virus, starting at time t. Q m (t) can be expressed as an integration of new infections 
during the period of missing doses (according to Eqn. 5): 

t+m-T d />\ 

njt)~K ave Jt)-S- / Kx)dx = 7(0 ■ ™; W • (e < W»>*»r _ i) (ii) 

The expected number of target cells that become infected by fully resistant mutant viruses, O m , 
is a product of the effective mutation rate from the partially resistant mutant to the fully 
resistant mutant (fi e ff) and the total number of cells infected by the partially resistant mutant 

(Q m ): O m (0 = A*^ -£2(0, as shown in Eqn.4. 
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Note that we track the population of newly infected cells to assess the risk of de novo generation 
of full resistance. This assumes implicitly that the fully resistant mutant is selected only when it 
enters a cell. This is a good assumption for DAAs that act on intracellular stages of the viral 
life-cycle, such as viral genome replication or assembly. However, in situations where the drug 
blocks viral entry into the cell, the mutant virus may have a selective advantage for entering a 
cell. Then the viral population should be tracked instead, but the results presented here still can 
be applied to drugs that block cell entry by multiplying with a simple scaling factor (Perelson and 
Nelson 1999). 

Stochastic-deterministic hybrid simulation of multiple strains of HCV 

We constructed a simulation model considering the dynamics of the baseline virus and all the 
potentially partially resistant mutants (see Supplementary Material). This simulation model 
follows a hybrid approach used previously to simulate the evolutionary dynamics of HIV(Ke and 
Lloyd-Smith 2012). It considers the dynamics of multiple strains of HCV deterministically (using 
ODEs) while treating the extinction and mutation processes as stochastic events (see 
Supplementary Material for detail). 

In the simulation, a patient is treated for a total period of 24 weeks. We generate two types of 
dosing patterns: random dosing and guided dosing. For the random dosing pattern, doses are 
missed in blocks of 1-3 days at times chosen randomly with equal probability during the 
treatment period. This probability is set as a constant in each run, but varied across runs such 
that different overall levels of adherence are generated. In each simulation, we assume that at 
least one-day treatment is taken immediately after each dose-skipping event (i.e. 1, 2 or 3 
consecutive missed doses), to ensure that two dose-skipping events do not occur consecutively 
(otherwise, longer blocks of doses would be missed than intended). For guided dosing, the 
procedure is the same as for random dosing, except that we ensure that: 1) doses are always 
taken during the high-risk window period predicted by our theory, and 2) after the window, a 
sufficient number of uninterrupted doses (calculated as N m ) are always taken immediately after 
missing doses, 3) if virus is not eradicated after the 24 weeks treatment period, a patient is 
treated with an uninterrupted number of doses as predicted by our theory, to ensure eradication 
of the virus. The outcome of the simulation at the end of the procedure is reported. 
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Figure 1. The impacts of suboptimal adherence on viral fitness. (A) A schematic illustrating 
how a non-preexisting mutant, A'B', fully resistant to a combination therapy involving two drugs, 
A and B, can be generated when adherence is suboptimal. Each black circle represents a mutant 
on the parameter space of resistance levels to A and B. AB, A'B and AB' are preexisting mutants 
that are non-resistant, resistant to A only and resistant to B only, respectively. Colored areas 
denote parameter regimes where mutants are fully resistant to the therapy (red), can grow if 
doses are missed (pink), and do not grow (blue). Note that the pink area can grow or shrink on 
the parameter space depending on the number of consecutively missed doses and drug PK/PD, 
and mutants lying in the pink area are 'partially resistant mutants'. (B) The dynamics of viral 
strains under treatment are determined by several factors: drug concentration, [D], which 
decreases with an increasing number of missed doses, m (upper panel); how viral replication is 
affected by drug [1-e; middle panel); and the relative number of target cells, /)(t) (lower panel). 
Upon effective treatment, h[t) increases to the infection-free level. (C) We integrate all these 
factors into a single fitness parameter, Reffif). Viral fitness increases as drug concentration drops 
(indicated by shades of green) and as target cell abundance rises (the blue arrow). Values of R e /f{t) 
can exceed 1, i.e. positive growth, if doses are missed after a period of effective treatment. 
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Figure 2. Suboptimal adherence prolongs treatment time needed to eliminate partially 
resistant mutants and increases the risk of de novo evolution of fully resistant strains. Two 

simulations assuming perfect adherence (dashed black lines) and imperfect adherence (solid red 
lines) are shown. In the simulation assuming imperfect adherence, single doses are missed at day 
5 and 15, and 3 consecutive doses are missed during days 22-24. (A) Drug concentration over 
time normalized by the maximum drug concentration C max - (B) The abundance of target 
hepatocytes, which rebounds after the initiation of combination therapy. (C) Viral fitness of the 
partially resistant mutant under consideration. Missing doses increases the value of R e ff, 
especially when multiple consecutive doses are missed or when h[t) has increased to high levels. 
(D) The dynamics of cells infected by the PM (on logio scale). The number of infected cells 
declines almost exponentially when doses are taken. Missed doses allow the number of infected 
cells to rebound. This means that an additional period of treatment is needed to suppress the 
mutant below the extinction threshold level, i.e. to achieve viral elimination. (E) The number of 
cells newly infected by the partially resistant mutants. Missed doses lead to substantial numbers 
of additional new infections, especially when 3 consecutive doses are missed. 
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Figure 3. There is a high-risk window early in treatment when missing doses is more likely 
to cause de novo resistance. (A) The changes in the risk of de novo resistance, O m , generated by 
a partially resistant mutant over time. The two sets of trajectories, A and B, differ in that the 
value of jUegrfor trajectory B is smaller by a factor of 10 5 (representing one additional nucleotide 
mutation) than the value set for trajectory A. Each set of trajectories shows the risk when the 
number of doses missed [m) is 1,2 or 3. (B) Dynamics of the two time-varying quantities in Eqn.4, 
i.e. the number of cells infected by the partially resistant mutant relative to the initial number 
before treatment (I(t)/I(0); blue dashed line), and the value of 0(t), green dotted lines, as 
shown in Eqn.4. Under effective treatment, the number of infected cells /(£j decreases 
exponentially, while the number of target cells rebounds to the infection-free level quickly, 
causing an increase in R aV e,m and thus 0(£). Together these changes cause <!>,„ to increase initially 
and then to decrease exponentially at longer times (as seen in panel A). 
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Figure 4. Adaptive treatment strategy improves treatment outcome substantially - a 
case when the risk of de novo resistance is low (with wild-type genotype- lb HCV at 
baseline). (A) Theoretical prediction of the treatment duration needed to eliminate each 
partially resistant mutant under perfect adherence (gray bar), and the maximum number of 
additional doses needed to compensate for missing doses, W m mox (colored bars). Green, blue 
and red denote results when 1, 2 and 3 consecutive daily doses are missed, respectively. The 
symbols next to the bars for /V m mox show the type of mutant investigated in panels (B,D,F). (B) 
Theoretical prediction of the risk of de novo resistance, O m , over time (as shown in Fig. 3d), 
for the three mutants with highest risks of generating fully resistant mutants. The dashed 
black line shows <J> m =0.01. (C) Treatment outcomes 1-3 days of doses are missed randomly. 
Colored areas denote the fractions of simulations with outcomes of viral relapse without full 
resistance (light blue) and viral clearance (gray). (D) Comparison between theory 
predictions and simulations of the number of cells infected by different mutants after 24 
weeks of treatment. (E) Treatment outcomes if adaptive treatment strategy is followed. The 
area above the black dashed line denotes the fraction of patients where virus is not cleared 
after 24 weeks' treatment. After 24 weeks, patients take the prescribed number of make-up 
doses without missing further doses. White areas denote adherence levels that are not 
allowed by the adaptive treatment strategy. (F) Same comparison as in panel (D) for the 
guided dosing simulation. 
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Figure 5. Adaptive treatment strategy prevents de novo resistance and improves treatment 
outcome substantially - a case when the risk of de novo resistance is high. Theoretical 
prediction and simulation for patients with the Y93H mutant virus (genotype-lb) at baseline 
under combination therapy of daclatasvir and asunaprevir. Thus, the mutants considered here all 
have the Y93H mutation. The theoretical predictions and simulation results are plotted in the 
same way as in Fig. 4. Dark red areas in panel (C,E) denote the fraction of patients with de novo 
full resistance to the combination therapy. Note that the fraction of patients with de novo 
resistance in the guided dosing scenario is very small (<0.1%). When doses are guided, so that 
mutant viral load does not rebound to the pre-treatment level, the theoretical prediction agrees 
well with simulation as shown in panel (F). 
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Supplementary Text 



1. HCV model and derivation of viral fitness (Ro) 

We first construct an ordinary differential equation (ODE) model to describe the long-term 
within-host dynamics of a single HCV strain under drug treatment. This model is based on an 
established model developed by Neumann et a/.(Neumann, Lam et al. 1998). It considers the 
dynamics of the target hepatocytes (H), the infected hepatocytes (I) and the HCV viruses (V). 
Guedj et al. have recently shown that the dynamics of viral loads during the first few days of 
treatment with direct acting antivirals (DAAs) are better described by a multi-scale model that 
considers intracellular dynamics of HCV RNAs(Guedj, Dahari et al. 2013). However, since here we 
are mostly interested in the longer-term dynamics of HCV infection, stretching for weeks or 
months, the simpler model by Neumann et al. is a good approximation. 



The ODE model describing the system is: 

dH 

—— = A — d-H — B-H-V 

dt 

dl 

- r P- H -v-s-, (S1) 

dV 

— ={l-e)- V -l-c-V 

Uninfected target hepatocytes are produced at constant rate, 1, and cleared at per capita rate d. 
The infection rate for uninfected hepatocytes is proportional to V, with infection rate constant 
p\The per capita death rate of infected hepatocytes is 5. The virions are cleared at per capita rate 
c. In the absence of drug treatment, viruses are produced from infected hepatocytes at rate p. 
Under drug treatment, the production of viruses from infected cells is reduced to rate {l-e)*p, 
where e is the efficacy of the drug of drug treatment. 



It has been shown that the number of target hepatocytes increases quickly after initiation of 
effective treatment (Rong, Dahari et al. 2010). In our model, this rebound rate is determined by 
the parameters 1 and d, and the number of target hepatocytes in the absence of HCV infection, 
H 0 , is determined by the ratio of this two parameters, i.e. H 0 =X/d. We chose values of 1 and d 
such that the target hepatocytes increase on a similar timescale to the results of Rong et o/.(Rong, 
Dahari et al. 2010), while keeping the total number of target hepatocytes constant in the 
absence of infection. To calculate R e /f{t) (in the main text), we set h{t) as h{t)=H{t)/H 0 , i.e. the 
normalized abundance of target cells. 

Based on Eqns. SI, we can calculate the reproductive number, R 0 , in the absence of drug, as: 

R 0 =H 0 -p-p/(8-c) (S2) 
We set /?o=10 in our main analysis, and this choice is in agreement with previous 
studies(Neumann, Lam et al. 1998; Rong, Dahari et al. 2010). 
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We assume there are 2*10 hepatocytes in an infected liver(Rong, Dahari et al. 2010). It has 
been shown that l%-50% of all hepatocytes are infected in chronically infected patient(Liang, 
Shilagard et al. 2009; Wieland, Makowska et al. 2014). Thus, we assume that only half of the 
total hepatocytes can potentially be infected in the absence of treatment, i.e. H 0 =l*10 11 , as in 
Rong et o/.(Rong, Dahari et al. 2010). The value of p is set such that the reproductive number R 0 
of the virus in the absence of drug is 10 (calculated above in Eqn. S2). The parameter values are 
listed in Table SI. 



Table SI. Parameter values in the HCV model. 



Parameters 


Values 


Unit 


References 


k 


1.95*105 


cells ml 1 day 1 


See text 


d 


0.15 


day- 1 


See text 




8.88*10-8 


ml day- 1 


Rong et a/. (Rong, 
Dahari et al. 2010) 


6 


0.15 


day- 1 


Neumann et 
al. (Neumann, Lam et 
al. 1998) 


P 


289.76 


day- 1 


See text 


c 


22.3 


day- 1 


Guedj et o/.(Guedj, 
Dahari et al. 2013) 



2. Combination therapy of daclatasvir and 
asunaprevir 



Pharmacokinetics/Pharmacodynamics 

In general, pharmacokinetics of HCV DAAs follow a characteristic pattern: drug concentration 
increases quickly to a peak level after dosing and then decreases exponentially until the next 
dose is administered. We use three pharmacokinetic parameters to describe this pattern: the 
time to reach peak concentration after taking the drug (r), the peak drug concentration {C max ) 
and the minimum drug concentration before the next treatment (C mm ). We assume the active 
drug concentrations in the liver are related to the drug concentration in the plasma (where data 
are measured) by a constant ratio r\. Then, the active tissue concentration of a drug, C(f), 
between dosing intervals can be described using the following equation: 

(C min + — t)-r, 0<t<r (s3) 

C max ■ exp (— w ■ (t — t)) -t] t < t < T 

where C mox and C min are the maximum and minimum drug concentrations in the plasma, 7" is the 

1 c • 

interval between two consecutive doses, and w = log^ 221 . The value of w is calculated 

T—T Cmax 

such that the drug concentration at the beginning is equal to the concentration at the end of a 
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single dose. 

The regimen used in clinical trials for the combination therapy of daclatasvir and asunaprevir is 
60mg once daily for daclatasvir and 200mg twice daily for asunaprevir(Lok, Gardiner et al. 2012; 
Pol, Ghalib et al. 2012). The value of liver-to-plasma ratio, rj, for daclatasvir is set as r\ — 0.094 
for daclatasvir as shown in recent work(Ke, Loverdo et al. 2014). The value of liver-to-plasma 
ratio for asunaprevir in human subject is still not clear, although this drug has been shown to 
have a large liver-to-plasma ratio in animal models(McPhee, Sheaffer et al. 2012). Clinical data 
from treated patients shows that mutant Q80L+D168V has EC 90 (the drug concentration suppress 
90% production) of 55 nM, and it is resistant to asunaprevir treatment in a patient (PT-29) with 
trough plasma concentration of 18-33nM(Karino, Toyota et al. 2013). This suggests that the 
tissue concentration of asunaprevir in the liver is at a similar level to its level in the plasma, and 
therefore, we have set r\ = 1.0 for asunaprevir. The parameters used for the pharmacokinetics 
of these two drugs are shown in Table S2. 



Table S2. Pharmocokinetic parameter values for daclatasvir and asunaprevir treatment used in 
the simulation model. 



Parameter 


Daclatasvir 
(Nettles, Gao 
etal. 2011) 
60mg QD 


Asunaprevir 

(Eley, 
Pasquinelli et 
al.2010) 
200mgBID 


Cmax 


1726 ng/ml 


268 ng/ml 


Cmin 


255 ng/ml 


35 ng/ml 


T 


1 day 


0.5 day 


X 


1.5 hour 


3.0 hour 




0.094 


1 



Since daclatasvir and asunaprevir act independently on NS5A and NS3 genes, here in the model, 
we calculate the inhibition of viral growth using Bliss independence(Bliss 1939): 

£;(t) = r^mi^m (S4) 

where C ooc (t) and C osu (t) are the active tissue concentrations of daclatasvir and asunaprevir, 
respectively, and EC 5 o_dac and EC S o_osu are the corresponding EC 50 values of for the viral strain 
under consideration. The average inhibition during the period when m doses are missed, £ a ve,m> 
is calculated numerically using Eqns. S3 and S4. Another model for pharmacological 
independence is Loewe independence(Loewe and Muischnek 1926). Changing the model to 
Loewe independence slightly changes our prediction about the fitness of each mutant under 
treatment, but does not alter the conclusion of the model. 



Characterizing preexisting mutants (PMs) and predicting the time needed to eradicate the PMs 

25 



Downloaded from http://biorxiv.org/ on September 18, 2014 



We first approximate the equilibrium level of cells infected by the baseline virus. The baseline 
virus is the viral strain that dominates the population before treatment, i.e. either the wild-type 
or the Y93H mutant in this study. Its equilibrium abundance before treatment can be derived 
from the single strain model in Eqn. SI, yielding: 

X c ■ d 

Under effective treatment, the population of infected cells declines exponentially, at a rate set 
by the half-life of the infected cell (l/d). Then, the time needed to eradicate the baseline virus 
under perfect adherence can be calculated as: 

terad.base = 0og/ 0 ~ log/ ext ) ■ j = (log(£ - j^) - log/ ext ) ■ | (S5) 

where l ext is the extinction threshold of infected cells below which the virus goes extinct. In our 
model, we set l ext =l/15000 copy/ml (assuming there are 15L of extracellular fluid(Rong, Dahari et 
al. 2010)). 

If the fitness of a mutant relative to the wild-type is l-s mut , then the frequency of resistant 
mutant virus before treatment can be approximated as(Ribeiro, Bonhoeffer et al. 1998): 

. _ Mmut . 
'mut — ~ '0 
$mut 

where fj, mut is the mutation rate from the baseline virus to the mutant. Then the time needed to 
eradicate the mutant virus, t era d jmu t, can be calculated as: 

terad.mut = (\ogI mut - log/ ext ) ■ | = (l<>g(^ " / 0 ) ~ !og hxt) ' \ ( S6 ) 

In the model, the mutation rates are set as 2.85*10" 5 and 1.5*10" 6 per infection cycle in infected 
cells for transitions and transversions, respectively (Loverdo et al, unpublished work). A mutant is 
considered as a preexisting mutant if the number of infected cells calculated for the mutant is 
above 1 copy in a patient, which corresponds to 1/15,000 copy/ml(Rong, Dahari et al. 2010). If a 
mutant strain has a higher replicative fitness than the wild-type virus (as measured in the 
replicon system), we set its relative fitness to 0.99 to ensure that the baseline virus is the 
dominant strain before treatment. This assumption is required by definition of the baseline 
strain, and reflects possible differences between fitnesses measured in replicon systems and in 
vivo. 

Characterizing fully resistant mutants 

Since daclatasvir and asunaprevir act independently on different target genes (NS5A and NS3, 
respectively), we define mutant viruses bearing resistance mutations, i.e. mutants show positive 
growth under therapy, against both daclatasvir and asunaprevir as potentially fully-resistant to 
the combination therapy. 

To characterize these mutants, we first find mutations that cause resistance, i.e. positive growth, 
under monotherapies of daclatasvir or asunaprevir, and we assume that daclatasvir and 
asunaprevir concentrations are the same with the corresponding concentrations in the 
combination therapy. For each mutant that have been reported to show higher resistance level 
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than the wild-type, we calculated the effective reproductive number under daclatasvir and 
asunaprevir monotherapies {R e ff,dac_max and R e ;iasu_max, respectively) when the target cell 
population is at its infection-free level (Table S3 and S4). We find that, if the wild-type virus is the 
baseline strain before treatment, the preexisting mutants that are resistant to daclatasvir 
monotherapy are L31M/V+Y93H, and the preexisting mutants resistant to asunaprevir 
monotherapy are D168A/V (Table S3). Then, combinations of these mutations are potentially 
fully resistant mutants, e.g. L31V+Y93H+D168V. If the Y93H mutant virus is the baseline strain 
before treatment, the preexisting mutants that are resistant to daclatasvir monotherapy are 
L31M/V+Y93H and L31V+Q54H+Y93H, and the preexisting mutants that are resistant to 
asunaprevir monotherapy is the same with the case when the wild-type virus is at baseline, i.e. 
D168A/V (Table S4). 



Table S3. The resistance profile of genotype lb mutants when the wild-type virus is at baseline. 



Mutantf 


EC50 to 


EC50 to 


Relative 


Reff,dac_max 


Reff,asu_max 




daclatasvir 


asunaprevir 


replication 








treatment* 


treatment* 


(1-Smut) * 






WT 


0.0026 


0.86 


1.0 


0.0004 


0.0716 


L31M 


0.0084 


0.86 


0.99** 


0.0011 


0.0708 


L31V 


0.0716 


0.86 


0.99** 


0.0096 


0.0708 


L31W 


0.2100 


0.86 


0.99** 


0.0281 


0.0708 


Q54H 


0.0032 


0.86 


0.83 


0.0004 


0.0594 


Y93H 


0.0621 


0.86 


0.27 


0.0023 


0.0193 


D168A 


0.0026 


109 


0.37 


0.0001 


1.6253 


D168V 


0.0026 


241 


0.29 


0.0001 


1.7968 


L31F+Y93H 


14.87 


0.86 


0.29 


0.4673 


0.0208 


L31M+Y93H 


18.23 


0.86 


0.70 


1.3248 


0.0501 


L31V+Y93H 


37.94 


0.86 


0.50 


1.5923 


0.0358 


Q54H+Y93H 


0.0243 


0.86 


0.22 


0.0007 


0.0157 


L31V+Q54H+Y93H 


48.74 


0.86 


0.99** 


3.6766 


0.0708 



t Bold names denote mutants that are preexisting in a patient. 

* Data taken from Fridell et al.(Fridell, Qiu et al. 2010) and McPhee et al.(McPhee, Friborg et al. 2012). 

** The fitness values of these mutants relative to the wild-type are set to 0.99, because the values measured in 

the replicon system were higher than 1 . 



Table S4. The resistance profile of genotype lb mutants when the Y93H virus is at baseline. 



Mutantf 


EC50 to 


EC50 to 


Relative 


Reff,dac_mctx 


Reff,asu_max 




daclatasvir 


asunaprevir 


replication 








treatment* 


treatment * 


(1-SmuO * 






Baseline-Y93H** 


0.0621 


0.86 


1.0 


0.0084 


0.0716 


L31F+Y93H 


14.87 


0.86 


0.29 


0.4673 


0.0208 


L31M+Y93H 


18.225 


0.86 


0.70 


1.3248 


0.0501 
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1 31 \/+YQ3H 


27 qac 
j i .j j j 


U. OD 


U. JU 


1 ^Q93 


U.UJJO 


Q54H+Y93H 


0.0243 


0.86 


0.22 


0.0007 


0.0157 


Y93H+D168A 


0.0621 


109 


0.37 


0.0031 


1.6253 


Y93H+D168V 


0.0621 


241 


0.29 


0.0024 


1.7968 


L31V+Q54H+Y93H 


48.735 


0.86 


0 99* ** 


3.6766 


0.0708 



t Bold names denote mutants that are preexisting in a patient. 

* Data taken from Fridell et a I .( Fridell, Qiu et al. 2010) and McPhee et al.(McPhee, Friborg et al. 2012). 

** We assume that the baseline virus, Y93H mutant, has the highest fitness, set its relative fitness as 1.0. 

*** The fitness value of this mutant relative to the wild-type is set to 0.99, because the value measured in the 

replicon system is higher than 1. 



3. Simulation of a hybrid multi-strain model 

We construct a simulation model considering the dynamics of the baseline virus and all the 
preexisting mutants that are shown in Table S3 and S4. This simulation model follows a hybrid 
approach used previously for simulating HIV evolutionary dynamics(Ke and Lloyd-Smith 2012). 
The model considers the dynamics of multiple strains of HCV deterministically, using ODEs, while 
treating the extinction and generation of mutants as stochastic events. The model tracks the 
population of uninfected and infected hepatocytes. Since the dynamics of viruses are much 
quicker than those of infected cells, we assume that the virus population is at quasi-equilibrium 
with respect to the dynamics of the infected cell population: the viral abundance is then given by 
V(t) — (1 — e) ■ p ■ I(t)/c. Then, the ODEs describing the dynamics of the multi-strain system 
become: 

dH B ■ p 

— = l-d-H- — -H-) (1 -£;)"/; 

at ct—i r . 

i=i (S7) 

dh B-p 

^ = ^. H . (1 _ £) ./ 8 . k 

dt c 

where H is the concentration of target hepatocytes, and /, is the concentration of hepatocytes 
infected by viral strain /'. 



The mutation process is treated stochastically. During the simulation, the ODEs are first 
simulated for a fixed time increment (At=0.01 day). At the end of each time increment, we 
approximate the number of cells newly infected by viruses from cells infected by the i th strain as 

— ■ H ■ (1 — £j) ■ /j ■ L ■ At, where L is the total volume of the liver. Of these newly infected cells, 

the number of cells in which the infecting viral lineage mutates from the i th strain to the j th strain 
can be drawn from a binomial distribution with probability of u^, which is the mutation rate 
from the i th to the j th strain. We then convert the number to concentration by dividing the 
number by L: 
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Mi < = BQ—!- -H- (1 -&)■/,■ L ■ M,u u )/L 
>j c j 

Each time step, the concentration of each strain is updated according to the values of A/y. We 

then check the number of cells infected by each strain. If the number is less than 1 copy per 

individual, we set it to 0 in the system, i.e. extinction. These procedures are iterated until all 

infected cells are extinct or simulation time exceeds 24 weeks (for random dosing pattern) or 

guided dosing period (for guided dosing). 

The dosing pattern is generated according to the procedure described in Online Methods. Once 
the dosing pattern is generated, the drug concentrations are calculated according to Eqn. S3. 



4. Sensitivity analysis 

In the analytical derivations, the parameters that determine the values of N m and O m are: the 
rate at which target hepatocytes become available under treatment (set by the values of 
parameters A and cf); the overall viral fitness, R 0 ; and the clearance rate of infected hepatocytes, 
8. We performed two rounds of sensitivity analysis, first testing how changes in these parameter 
values impact the clinical outcomes predicted by our theory, and then testing the robustness of 
our adaptive treatment strategy to changes in these parameter values. 

Sensitivity of predicted clinical outcomes to changes in 
parameter values 



The overall viral fitness, R 0 

The viral fitness parameter, R 0 , has several impacts on the predictions of the model. To evaluate 
the impact of changes in R 0 on the time needed to eradicate the virus, we can substitute the 
expression of R 0 into Eqn. S5 as: 

terad.base = (log(^ " (l ~ ^)) ~ log I ext ) ■ | (S8) 

From this equation it can be seen that, if A and 8 are kept constant, higher R 0 leads to a higher 
level of infected hepatocytes before treatment, and thus, a longer time needed to eradicate the 
virus. However, this increase may be small, because t em d,base changes in proportion to the 

i 

logarithm of 1 . 

#0 

The value of N m scales linearly with R 0 (Eqn. 2 in the main text). Thus, the higher the value of R 0 , 
the more doses needed to compensate for any missed doses. 

The value of O m increases almost exponentially with an increase in R 0 (Eqn. 3 in the main text). 
Therefore the risk of generating fully-resistant mutants increases drastically as R 0 rises (see Fig. 
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Sla,b). This is because the viral population of a mutant increases exponentially if R e jf,m becomes 
greater than 1, which leads to an exponential increase in the risk of generating fully resistant 
mutants. 

The clearance rate of infected hepatocytes, 8 

The clearance rate of infected hepatocytes, 8, influences our theoretical predictions in two ways. 
First, the time needed to eradicate the virus depends linearly on the inverse of the clearance rate, 
1/8 (Eqn. S5 and S6). The quicker the clearance rate of infected hepatoctyes, the shorter time 
needed to eradicate the virus. Second, changes in the value of 8 affect the risk of generating fully 
resistant mutants when doses are missed (Eqn. 3 in the main text). If we let the fitness of a virus, 
/?o, unaffected by changes of 8, larger values of 8 lead to shorter half-lives of infected 
hepatocytes, and higher risk of generating fully-resistant mutants, because the viral lineage 
undergoes more rounds of replication during a fixed dosing period (Fig. Sic). 

The rate at which target hepatocytes become available under treatment 

N m is linearly dependent on the level of target hepatocytes (Eqn. 2 in the main text), and thus, 
slower rebound of target hepatocytes would decrease the number of compensatory doses 
needed if doses are missed before the target hepatocyte population rebounds back to its 
infection-free level. 

The rebound rate of target hepatocytes also impacts O m , through its influence on the effective 
viral fitness, R e ff,m (Eqn. 3 in the main text). In our model, the rebound rate is set by the values of 
parameters A and d. If we keep the number of hepatocytes before treatment constant, by 
keeping the ratio of A over d constant, then a slower rebound rate (lower value of both A and d) 
results in a substantially reduced rate of generating de novo resistance (Fig. Sid). 

Robustness of adaptive treatment strategy to variations in 
parameter values 

We first tested the robustness of adaptive treatment strategy to lower or higher values of R 0 
(/? 0 =5 and /?o=15, as opposed to /? 0 =10 for our main results). When /?o=5, both the number of 
compensatory doses (N m ) and the potential to generate de novo resistance (<l> m ) decrease 
substantially, as predicted by our model (Fig. S6-S7). This leads to a shorter high-risk window 
period when mutant Y93H is the baseline strain (Fig. S7) and lower adherence levels are required 
to eradicate the virus. When /?o=15, we observe the opposite pattern: higher adherence levels 
are required to eradicate the virus and there is a longer high-risk window period when mutant 
Y93H is the baseline strain (Fig. S8-S9). 

We then tested how well our adaptive treatment strategy works when the half-life of infected 
hepatocytes is shorter. As shown in Fig. S10 and Sll, the time needed to eradicate the virus 
decreases substantially, to 12 weeks of effective treatment. Our adaptive treatment strategy 
improves clinical outcome especially when Y93H mutant is at the baseline (Fig. Sll). 
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In general, we find our adaptive treatment strategy is robust against variations in key parameter 
values. Under all parameter values tested, the adaptive treatment strategy delivers 
substantially better patient outcomes than random dosing with the same overall adherence 
levels. 
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Supplementary Figures 
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Figure SI. Sensitivity analysis of the risk of de novo resistance to variations of key parameter values, 
Ro.mut (panels A,B), 8 (panel C) and the rate of recovery of target cells upon treatment (panel D). In 

each panel, the trajectories show how the risk of de novo resistance (Log^^m) changes over time if 
adherence is perfect. Figures are plotted using the same parameter settings as trajectories 'a' in Fig. 3 
in the main text, except that Ro,mut=5 in panel A, R 0 , mu t=15 in panel B, <5=0.5 in panel C and 
a=1.95*10 4 , d=0.015 in panel D. In the main results, i.e. Fig. 3, the parameter values used are 
Ro, mu t=10, <5=0.15, A=1.95*10 5 , d=0.15. 
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wild-type at baseline 




50 100 
Days under treatment 

Figure S2. The predicted number of additional days of doses needed to compensate for the 
first instance of missing 1 (green lines), 2 (blue lines) or 3 (red lines) consecutive days of doses 
(maximum N m for all partially resistant mutants), and the high-risk window period of de novo 
resistance (shaded area; O m >0.01 for any of the partially resistant mutants as shown in Fig. 4B). 

(A) Predictions for patients with the wild-type virus at baseline before treatment. The areas 
below the curves are white, indicating that the risk of de novo resistance is always low. (B) 
Predictions for patients with the Y93H mutant virus at baseline. The initial increases of the 
number of days of compensating doses are due to the increase of the number of target cells 
upon treatment, and the sudden drops during later periods of treatment are due to the 
elimination of particular partially resistant mutant lineages. Note that these curves are 
calculated under the assumption that adherence is perfect except for the 1-3 days of missed 
doses being considered, i.e. it is a prediction for the first instance of missed doses. For cases 
where multiple instances of missed doses have occurred, one needs to calculate the values of N m 
and <J> m for each mutant based on the adherence pattern, and then integrate them together by 
choosing the highest values of N m and <5 m for those mutants. 
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Theory prediction Theory prediction 

Figure S3. Theory correctly predicts the number of cells infected by two mutants (L31V and 
L31M+Y93H) generated in the hybrid model simulation when doses are missed randomly (panels 
A,B) or guided by adaptive treatment theory (panels C,D) for patients with the wild-type virus at 
baseline. L31V and L31M+Y93H are the two most likely mutants that generate full resistance. The 
axes are the theory prediction (x-axis) and model simulation (y-axis) of the Log 10 of the number of 
mutants, which are calculated as the cumulative numbers of Log 10 ^ m (t)lji mut for all missed doses. 
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Figure S4. The theory correctly predicts the number of cells infected by mutant viruses at the end of 
24-weeks' treatment in the hybrid model simulation when adherence is greater than 70% (vertical 
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dashed lines) and doses are missed randomly, for patients with the Y93H mutant virus at baseline. 

The y-axis shows the log 10 difference between the theory prediction and the model simulation at the 
end of the 24-weeks' treatment. Note that when adherence is lower than 70%, the population of 
infected cells grows to high levels close to the pre-treatment level, where further growth is curtailed 
by target cell limitation. As a result, the theoretical prediction overestimates the number of cells 
infected by the mutant virus significantly because we assume the number of target cells is not limited. 
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Figure S5. Comparison between theory prediction and simulation of the numbers cells infected by 
three mutants (L31M+Y93H, L31V+Y93H, L31V+Q54H+Y93H ) generated when doses are missed 
randomly, for patients with the Y93H mutant virus at baseline. (A,B,C) The axes are the theory 
prediction (x-axis) and model simulation (y-axis) of the Log 10 of the number of cells infected by 
different mutants, which are calculated as the cumulative numbers of Log 10 O m /,u m ut for all missed 
doses (L31M+Y93H in panel A; L31V+Y93H in panel B; L31V+Q54H+Y93H in panel C). (D,E,F) Our 
theory prediction is accurate for adherence greater than 70%, but overestimates the number of cells 
infected by the mutant virus significantly when adherence is lower than 70%, for the same reason as 
explained in the legend of Fig. S4. 
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Figure S6. Theory correctly predicts the number of cells infected by three mutants (L31M+Y93H, 
L31V+Y93H, L31V+Q54H+Y93H) generated in the hybrid model simulation when doses are guided by 
adaptive treatment theory, for patients with the Y93H virus at baseline. (A,B,C) The axes are the 
theory prediction (x-axis) and model simulation (y-axis) of the Log 10 of the number of mutants 
(L31M+Y93H in panel A; L31V+Y93H in panel B; L31V+Q54H+Y93H in panel C), which are calculated as 
the cumulative numbers of Log 10 O m (t)//^ mut for all missed doses. (D,E,F) the Log 10 differences 
between theory prediction and model simulation as shown in panels (A,B,C). Note that our theory 
agrees very well for mutants L31M+Y93H and L31V+Y93H. For mutant L31V+Q54H+Y93H, the 
stochastic extinction and appearance of this mutant generates stochastic deviations of the simulation 
from theory prediction. 
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Figure S7. The impact of lower viral fitness (R 0 =5) on treatment outcomes and adaptive 
treatment strategies of combination therapy with daclatasvir and asunaprevir, with the 
wild-type virus at baseline. Panels A-F show the same plots as Fig. 4 in the main text, except that 
the fitness parameter R 0 for the wild-type virus is assumed to be 5. The treatment outcome 
improves for all scenarios for this lower viral fitness (compare with Fig. 4 in the main text). Using 
the adaptive treatment strategy prevents viral relapse and de novo resistance if overall 
adherence is greater than 60% (panel E). Panel F is empty because all patients are cleared of 
infection after 24 weeks. 
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Figure S8. The impact of lower viral fitness (R 0 =5) on treatment outcomes and adaptive 
treatment strategies of combination therapy with daclatasvir and asunaprevir, with the Y93H 
virus at baseline. Panels A-F show the same plots as Fig. 5 in the main text, except that the 
fitness parameter R 0 for the wild-type virus is assumed to be 5. The treatment outcome improves 
for all scenarios for this lower viral fitness (compare with Fig. 5 in the main text). Using adaptive 
treatment strategy reduced the risk of de novo resistance (panels E). Our theory correctly 
predicts the number of infected cells in a patient at the end of 24 weeks' treatment. 
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Figure S9. The impact of higher viral fitness (R 0 =15) on treatment outcomes and adaptive 
treatment strategies of combination therapy with daclatasvir and asunaprevir, with the 
wild-type virus at baseline. Panels A-F show the same plots as Fig. 4 in the main text, except that 
the fitness parameter R 0 for the wild-type virus is assumed to be 15. The treatment outcome 
improves for all scenarios for this lower viral fitness (compare with Fig. 4 in the main text). Our 
theory correctly predicts the number of infected cells in a patient at the end of 24 weeks' 
treatment. 
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Figure S10. The impact of higher viral fitness (R 0 =15) on the treatment outcomes and adaptive 
treatment strategies of combination therapy of daclatasvir and asunaprevir with the Y93H 
virus at baseline. Panels A-F show the same plots with Fig. 5 except that in the analytical 
derivation and model simulation, the fitness parameter R 0 for the Y93H mutant virus is assumed 
to be 15. The risks of viral relapse and de novo resistance become higher when the viral fitness, 
R 0 , is higher. Using adaptive treatment strategy can prevent de novo resistance and improve 
treatment outcomes (panels E and F). Our theory correctly predicts the number of infected cells 
in a patient at the end of 24 weeks' treatment when doses are guided (panels F). Our theory 
does not predict the number of infected cells at the end of treatment well, when doses are 
missed randomly and the adherence is low. This is because, when adherence is low, the viral load 
often rebounds back to the pre-treatment level, where it is limited by target cell availability. 
This phenomenon is not included in our theory, which overestimates the number of viruses as a 
result. 
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Figure Sll. The impact of higher viral clearance rate (5=0.5) on the treatment outcomes and 
adaptive treatment strategies of combination therapy of daclatasvir and asunaprevir with the 
wild-type genotype lb virus at baseline. Panels A-F show the same plots with Fig. 4 except that 
in the analytical derivation and model simulation, the viral clearance rate, 6, is assumed to be 0.5 
instead of 0.15 in Fig. 4 (but note that R 0 for the viruses is kept the same). When the viral 
clearance rate increases, it takes less time to eradicate the virus from a patient. However, when 
doses are missed, the population of mutant viruses expands more quickly, because the half-life 
of the infected cells is shorter and thus it undergoes a higher number of replication generations 
during the period of missed doses. Using the adaptive treatment strategy can prevent viral 
relapse and de novo resistance and improve treatment outcome (panels E). Our theory correctly 
predicts the number of infected cells in a patient at the end of 24 weeks' treatment when doses 
are guided (panel D). 



43 



Downloaded from http://biorxiv.org/ on September 18, 2014 




Adherence No. of infected cells in Log [<} (Theory prediction) 

Figure S12. The impact of higher viral clearance rate [3=0.5) on the treatment outcomes and 
adaptive treatment strategies of combination therapy of daclatasvir and asunaprevir with the 
Y93H virus at baseline. Panels A-F show the same plots with Fig. 5 except that in the analytical 
derivation and model simulation, we assume the viral clearance rate, d, is 0.5 instead of 0.15 (but 
note that R 0 for the viruses is kept the same). As seen in Fig. Sll for the scenario with the 
wild-type virus at baseline, it takes less time to eradicate the virus from a patient for this higher 
viral clearance rate. However, when doses are missed, the population of mutant viruses expands 
more quickly, increasing the risk of viral relapse and de novo resistance. Using adaptive 
treatment strategy can prevent viral relapse, de novo resistance and improve treatment outcome 
(panels E and F). Our theory does not predict the number of infected cells at the end of 
treatment well, when doses are missed randomly and adherence is low. This is because, during 
the time period when doses are missed, the rebound of the viruses is quicker when 6 is higher 
(because the viral generation time is shorter). When adherence is low, the viral load often 
rebounds back to the pre-treatment level, where it is limited by target cell availability. This 
phenomenon is not included in our theory, which overestimates the number of viruses as a 
result. 
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